function S = NeweyWest(hFunc,qLag)
% The required correlation matrices
T      = size(hFunc,1); %Number of observations
numMom = size(hFunc,2); %Number of moments
GAMA_array = zeros(numMom,numMom,qLag);
GAMA0 = CorrMatrix(hFunc,T,numMom,0);
if qLag > 0
    for i=1:qLag
        GAMA_array(:,:,i) = CorrMatrix(hFunc,T,numMom,i);
    end
end

% The estimate of S
S = GAMA0;
if qLag > 0
    for i=1:qLag
        S = S + (1-i/(qLag+1))*(GAMA_array(:,:,i) + GAMA_array(:,:,i)');
    end
end
end
